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Abstract 

We present a finite element method for the Stokes equations involving two immiscible incom- 
pressible fluids with different viscosities and with surface tension. The interface separating the 
two fluids does not need to align with the mesh. We propose a Nitsche formulation which al- 
lows for discontinuities along the interface with optimal a priori error estimates. A stabilization 
procedure is included which ensures that the method produces a well conditioned stiffness matrix 
independent of the location of the interface. 

Keywords: cut finite element method, Nitsche, Two-phase flow, Stokes, discontinuous viscosity, 
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1. Introduction 

A large number of real world phenomena exhibit strong or weak discontinuities. The application 
we have in mind is multiphase flow with kinks in the velocity field and jumps in the pressure field 
as well as in physical parameters, such as viscosity, across interfaces that evolve with time and 
may undergo topological changes. When simulating such phenomena, discontinuities can occur 
anywhere relative to a fixed background mesh. Unfortunately, standard finite element methods, 
as well as finite difference schemes, do not accurately model discontinuities that are not a priori 
fitted to the mesh. However, letting the mesh conform to the interfaces requires remeshing as the 
interfaces evolve with time, and leads to significant complications when topological changes such 
as drop-breakup or coalescence occur. Numerous strategies have been proposed to handle these 
difficulties. 

A common strategy has been to regularize the discontinuities [T|. However, this strategy has 
the drawback that it gives reduced accuracy near the interfaces and consequently requires a very 
fine mesh in these regions. Methods that allow for discontinuities along interfaces that do not align 
with the mesh, and hence avoid both regularization and remeshing processes, have become highly 
attractive and significant efforts have been directed to their development, see e.g. [3l [4j [5J [6] . In 
the finite element framework the extended finite element method (XFEM) , where the finite element 
space is enriched so that discontinuities can be captured [7], has become a popular alternative. 
However, in XFEM the conditioning of the problem is sensitive to the position of the interface. 
Whenever the interface cuts an element in such a way that the ratio between the areas/volumes on 
one and the other side of the interface becomes very large, the system may become ill-conditioned. 
In such cases, iterative linear solvers may breakdown. For unsteady problems it is not unusual 
that such situations occur. In [8], this problem is addressed by neglecting basis functions in the 
XFEM space that have very small support and may cause ill-conditioning. The criterion for the 
selection of basis functions to neglect then has to be chosen carefully so that accuracy is not lost. 

An alternative to the XFEM approach is based on an extension of Nitsche's method [9] for 
the weak enforcement of essential boundary conditions. This approach was first proposed for an 
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elliptic interface problem in [10] and later for a Stokes interface problem in [3j. The idea is to 
construct the discrete solution from separate solutions defined on each subdomain separated by 
the interface and at the interface enforce the jump conditions weakly using a variant of Nitsche's 
method. By choosing the coefficients in the Nitsche numerical fluxes locally on each element and 
letting them depend on the relative area/volume of each side of the interface the unfitted finite 
element methods in |10[ |3j can allow for discontinuities internal to the elements with optimal 
convergence order. However, these methods suffer from ill-conditioning just as XFEM. In [TT] and 
later in [T2] a stabilization of the classical Nitsche's method for the imposition of Dirichlet boundary 
conditions on a boundary not fitted to the mesh was considered for the Poisson problem and for 
the Stokes equation, respectively. In these methods the stabilization is applied in the boundary 
region and optimal convergence order and well conditioned system matrices are ensured. For the 
elliptic interface problem other stabilization strategies have been suggested as remedies to the 
ill-conditioning problem, see e.g. [131 [Uj. We also refer to [T5] for implementation aspects in three 
dimensions and |16| for extensions to higher order elements. 

In this paper, we propose a robust finite element method which offers a way to accurately 
and with control of the condition number of the system matrix solve the Stokes interface problem 
involving two immiscible fluids with different viscosities and surface tension. The model consists 
of the incompressible Stokes equations in two subdomains each occupied by a fluid. Differences in 
viscosity between the fluids and the surface tension force poses jump conditions at the interface 
separating the fluids. The finite element method we propose for this Stokes interface problem 
enforces the jump conditions at the interface weakly with weighted coefficients in the Nitsche 
numerical fluxes. We also suggest slight changes to the variational formulation in [3] to reduce 
spurious velocity oscillations and we include stabilization terms both for the velocity and the 
pressure that guarantee a well conditioned system matrix. The stabilization terms are consistent 
least squares terms controlling the jump in the normal gradient across faces between elements in 
a neighbourhood of the interface. Using the stabilization terms we prove an inf-sup condition and 
that the resulting stiffness matrix has optimal conditioning. We also prove the inf-sup stability 
of the method under the condition that the mean value of the pressure in the entire domain is 
fixed, in contrast to the inf-sup result in [3J which is based on the more restrictive condition that 
the mean values in each of the two subdomains are fixed. Our inf-sup result is also uniform with 
respect to the jump in the viscosity. The proposed method is simple to implement as it uses 
standard continuous linear basis functions with changes only in the variational form. 

The outline of this paper is as follows. In Section 2 we formulate the Stokes system and the 
finite element method. In Section 3 we prove that the method is of optimal convergence order. 
In Section 4 we prove that the condition number is (D(h~ 2 ) independent of the position of the 
interface relative to the mesh. Finally, in Section 6, we show numerical examples in two space 
dimensions and compare the method to existing methods. We summarize our results in Section 7. 

2. The interface problem and the finite element method 

We consider a problem consisting of two immiscible fluids separated by an interface, with the 
flow described by the incompressible Stokes equations. The Stokes system is a standard model for 
creeping viscous flow. In this section we present the equations and a finite element method for 
their approximate solution. 

2.1. The two- fluid incompressible Stokes equations 

Let O be a bounded domain in R 2 , with convex polygonal boundary <9f2. We assume that two 
immiscible incompressible fluids occupy subdomains Oj C £1, i = 1,2 such that fl = tti U SI 2 and 
fii =0 and that a smooth interface defined by T — dfii Hdf^ separates the immiscible fluids. 

We consider the following Stokes interface boundary value problem modeling two fluids with 
different viscosity and with surface tension: find the velocity u : fl — > M 2 and the pressure 
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p : fl -> E such that 

-V • (/xe(u) - pi) = / mfiiUn 2 , (2.1a) 

V-m = onOiUfi 2! (2.1b) 

[u] =0 on T, (2.1c) 

[n • (Ate(u) - pi)] = ctkti on T, (2. Id) 

u = g on <90. (2.1e) 

Here e(u) = (Vtt + (Vtt) T )/2 is the strain tensor, = 2/ij > on f2j,« = 1,2 is a piecewise 
constant viscosity function on the partition Qx U 0,2, f G [L 2 (f2)] 2 and g G [i/ 1 / 2 (9il)] 2 are given 
functions, cr is the surface tension coefficient, k is the curvature of the interface, n is the unit 
normal to T, outward-directed with respect to fii, and [a] = (&i — a2)|r, where dj = ala^i = 1,2. 

We assume a constant surface tension coefficient a. Hence, from equation (2. Id I, it follows 
that across the interface the shear stress is continuous, i.e., 

\n ■ fJ,e(u)j ■ t = 0. (2.2) 

Hereinafter, we also denote the outward directed unit normal vector on dfl by n. We assume 
global conservation of mass, that is 

/ g ■ nds = 0. (2.3) 

JdQ. 

We will use the notation (■,-) u f° r the L 2 (lu) inner product on lo (and similarly for inner 
products in [L 2 (uj)] 2 and [L 2 (w)] 2x2 ). We let |H| S , W and \v\ S}U denote the Sobolev norms and 
seminorms associated with the spaces H s (ui), respectively. For the convergence analysis we will 
assume that u G {H 2 (Sli U 2 )] 2 — W and the pressure space is 

V = {pe H 1 ^ U fia) : (m"V, lKun, = 0}. (2.4) 

2.2. Mesh and assumptions 

Let Kh be a triangulation of f2, generated independently of the location of the interface T. 
Introduce the set of all element faces JF associated with the mesh ICh, the set of all elements that 
intersect the interface 

/C r = {K G IC h : \T f)K\ > 0}, (2.5) 
and the set of all elements on the boundary 

}Cau = {K e fC h : \dQ nK\> 0}. (2.6) 

Define meshes on the subdomains f2j,i = 1, 2, as follows 

K h ,i = {K G IC h : 3 x, y, x £ y G K n H,} (2.7) 

and let 

Q M = (J K, uh ti = (J K, i = l,2. (2.8) 

Ke/Ch,i KeK h ,i,Kcn z 

Note that the interface T is allowed to intersect the elements in ICh and /C/jj and that elements 
intersected by the interface are both in Q^j and Qh,2- 

Let /Cr be the set of all elements that share two of its faces with elements in /Cr, let 

J- r<i = {F G T : F c <9AT, AT G /C r U £ r , F n fi, + 0, F n dfl = 0} (2.9) 

be the set of faces of elements in /Cr U /Cr, that have a nonempty intersection with fij, and are 
not on the boundary dfl, and let 

Uh,i = &h,i\ [j K ' i = 1 ' 2 > ( 2 - 10 ) 
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Figure 1: A triangulation of the domain Q. The interface T is a circle separating the subdomains Qi and Q2- 
Left panel: Oj, i (shaded triangles), 2 (non shaded triangles), and edges in Jpi (thick lines). Right panel: fl^ 2 
(shaded triangles), ujj (non shaded triangles), and edges in Jr,2 (thick lines). 



contain the elements in Qh,% that are not in ]Cr- We modify the set u>h,i m order to ensure that 
there are no elements in uihj that has an element with two edges on the boundary dbJh,i which 
is a necessary condition for the inf-sup condition condition to hold on LJh,i> see Lemma 3.10 See 
Fig. [Tjfor an illustration of the sets £lh,i-> u h,ii an d J~r.ii i=l,2 in a two-dimensional case. 
We make the following assumptions: 

• Assumption 1: We assume that the triangulation is quasi-uniform, i.e., there exists con- 
stants ci and C2 such that 

dh <h K < c 2 h \/k e lC h , (2.11) 

where Kk is the diameter of K and h = max^ /i^-, and that there are no elements with two 
edges on the boundary d£l. 

• Assumption 2: We assume that T either intersects the boundary dK of an element K £ ICr 
exactly twice and each (open) edge at most once, or that TDK coincides with an edge of 
the element. 

• Assumption 3: Let Tk,h be the straight line segment connecting the points of intersection 
between T and dK. We assume that Yk — T n K is a function of length on r^ ^; in local 
coordinates: 

?K,h = {(£, fj) ■ < £ < \T k , h \,r) = 0} (2.12) 

and 

r K = {(f, v) ■ < C < \r kth \,v = (2.13) 

• Assumption 4: We assume that for each K € ICr there are elements K % C $1^, i = 1, 2 such 
that Kn K 1 ^ 0. 

• Assumption 5: We assume that the mesh coincides with the outer boundary dft. 

Assumptions 2-3 essentially state that the interface is well resolved by the mesh. Except that 
we in the second assumption also allow the interface to be aligned with a mesh line these two 
assumptions are as in [10 . Assumption 4 states that each K £ ICr shares a face or at least a 
vertex with an element K 1 C f^i and an element K 2 C f^- For each Qi this is the same assumption 
as in IT2"1. 
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2. 3. Penalty parameters and averaging operators 

In our finite element method we will make use of the following averaging operators 

{a} = Kxdi + K 2 a 2 , (a) = n 2 ai + K\a 2 , (2-14) 

where the weights K\ and k 2 are real numbers satisfying k\ + k 2 = 1. We consider the following 
two cases, in accordance with the intersection options of T: 

• For each element A G /Cr such that T intersects the element boundary exactly twice, we 
have \K D f2,| = a i K h 2 K and |T n K \ = ^xhx for some cti t K, Ik > 0, and we define 

«1X = , > K 2\K = : , (2.15) 

a penalty parameter 

Ark = ^, (2.16) 
fix 

where 

m =D{n}+ (l + BhKHU* AG (0,1), B,^>0, (2.17) 
and if A is also an element in ICgn we define a penalty parameter 

A 8 dxn,, i ^( Gft+ 2(1 + F j WJ ' ). .Efo, ^ V (2.18, 
h K \ La K J \ l + max.\N K i)J 

Here, F,G > 0, A" 1 is the closest element to A that is completely in fij, and N K , is the 
number of elements in ICqq PI /Cr that are not completely in Qi and have the same element 
K l as the closest element completely in f2j. Furthermore, c q = max (j^ry, |fT/i)' wriere F 
are the faces that have to be crossed to pass from A to K l and \d£l K\ ~ "fdn^hx- 

• If r n A coincides with an element edge, then TDK will also be an edge of another triangle 
T G fCr and TDK = TOT. We may without loss of generality assume that J 1 C Sli 
and A c fi 2 . We may write \T\ = a T h\ and \K\ = a> K h K for some ar, > 0, and 
|T n A" I = |r PI T\ = ^k^k = "frhr for some jk, It > 0. We define 

Klk= 72— , «2k = " : 7—^2. (2.19) 



a penalty parameter 

Ar I u = 

|r n a 

where 



Ark = (2-20) 



m = D lK { f ,}+— (1+^)^2 AG (0,1); B,D>0, (2.21) 

A ((piaK)/T K + (a i 2 Q! t)/7t) 

and if A is also an element in /Can we define the penalty parameter 

Aaokno, = ^^+ (1 + f J^^ V Se( V , ^4 c V (2-22) 



Kk V Aa/f / \ l + max(A^^i) 

Here F, G > 0. For all the elements A G /Can that are not intersected by the interface we choose 
Aanknfii as in (2.22). Note that < m < 1 and K\ + k 2 = 1 in both cases. 
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Figure 2: Illustration of the domains, meshes, and spaces in a one-dimensional model case. 

2.4- The finite element method 

Let K-h be a triangulation of f2 that satisfies all the assumptions in Section 2.2 We let V/i )jU 



be the space of continuous piecewise linear polynomials defined on K-h with (fj, 1 qf l , l)niun 2 = 
Vqh & Vh,fj, an d we let 

V h =V M eV M (2.23) 

be our pressure space where 

V h ,i=V h>li \n hii , i = l,2, (2.24) 

i.e. the spaces of restrictions to tth,i and fth.2 of functions in Vh,n- A pi £ Vj consist of a pair 
(Ph,i,Ph,2), with ph,i e Vh,i- For an illustration in a one-dimensional model case, see Fig. [2} Note 
that ph G Vh is double valued on elements in JCr and is allowed to be discontinuous at the interface 
T. We define (p h , g , h)n 1 un 2 = Yn=i(Ph,i,qh,i)ni- 

In the same way as above we construct the velocity space but on a uniform refinement of Kh, 



denoted by K-h/2- The triangulation AC/,/2 also satisfies all the assumptions in Section 2.2 and we 



define all the parameters in Section 2.3 on the mesh To construct the velocity space we 

let Wh,o be the space of vector valued continuous piecewise linear polynomials on /C/j/2- We will 
impose the Dirichlet conditions weakly. Thus, there are no special boundary restrictions at dfl on 
the velocity space. We define 

W M = W M |fi M , * = 1,2 (2.25) 

and we let 

Wh = W h ,i®W h , 2 . (2.26) 

Each Uh € Wh consists of a pair (1*^,1, 14^,2)1 where Uh,i € Wh,i- 

We drop the subscript h and h/2 and denote both the velocity and the pressure mesh by K,. 
An element K G K, is an element in whenever we have terms involving functions in the 

velocity space Wh and it is an element in K-h whenever we have terms involving only functions in 
the pressure space. 

We propose the following Nitsche method: find (uh,Ph) & Wh x Vh such that 
A h (u h ,p h ;v h ,qh) + e u J(u h .v h ) +e p J{ph,qh) = L h (v h ), ^(v h ,q h ) e W h x V h - (2.27) 
Here Ah(-] •) is a bilinear form defined by 

A h (w,r;v,q) = a h (w , v) + b h (w , q) - b h (v,r), (2.28) 



G 



where 



a h (w,v) = (^e(w),e(v))n lU n 2 

- ({n ■ fxe(w)}, M) r - (M , {n ■ »e(v)}) r 

- (n ■ fj,e(w),v) d n -(w,n- ^,e(v))gn 
+ A r (M , W)r + AanO, v) d n 

b h {w, q) = (V • w, q)n lU n 2 - (In ■ wj , {q}) r -(n-w, q) dn 
b h (w, q) = -(w, Vq)n 1 un 2 + ([?] , (w • n)) r , 

and Lft(-) is a linear functional defined by 

L h (v) = (/, v)n + {an, (v ■ n))r + (g, XqV ~ n ■ fie(v) + nq) dn . 



(2.29) 
(2.30) 
(2.31) 

(2.32) 



In equation (2.27) e u and e p are positive constants and the stabilization terms are defined as 



J(j>h,qh) 



-U3/ 



»=1 FeFr.i 



[n F ■ VphA F , {riF ■ Vg ft 



(2.33) 



and the component wise extension for vector valued functions Uh,i = (u\ i5 4 ) 



J{u h ,v h 



2 2 

EE E 

J= l i=i FeJr,, 



Mi/l S ( 



v< 4 



JlJ? 



)f- 



(2.34) 



We choose s in different ways depending on if the interface cuts the domain boundary or not. We 
take s = 1 only for those faces that has to be crossed to pass from an element on the boundary 
K G tCgn, K (f. Q, , K n i ^ 0, to the closest element K l c Qj, otherwise s — 3. We have 
employed the following notation for the jump in a function v at an interior face F 



(2.35) 



where v ± (x) — lim t _j. + v ( x T trip), for x 6 F, and nj? is a fixed unit normal to F. 



Remark. By integrating the bilinear form 6/j in equation ( |2.30 ) by parts we get bh in equa- 
tion (2.31). We recommend to use bh in equation (2.31) in simulations since it results in reduced 



spurious velocities, see the remark in Example 2 of Section [5j 

Remark. The computation of the bilinear forms ah and bh require integration over Qi C &h,i- 
Thus, for elements cut by the interface T the integration should be performed only over parts 
of the elements. The functions Ph % and Uh % are defined on the larger subdomains £lh,i a- n d the 



stabilization terms (2.33) and (2.34) ensure well defined extensions from f2, to ft 



h , i • 



Remark. If the interface does not cut the boundary of the domain f2, the stabilization term 
J(ph,qh) in equation (2.33) is sufficient to prove the inf-sup stability of the method. However, to 



have control of the condition number of the system matrix independently of the position of the 
interface relative to the mesh both the stabilization terms J(ph,qh) and J(uh,Vh) are needed, 
see Fig. [9j The sets J-r.i, i=l,2 in equation (2.33) are defined for the pressure mesh TCh and the 



two sets J-r,i, i=l,2 in equation (|2.34|) are defined for the velocity mesh an d are different 



The face F is always a full face in the underlying mesh. Similar stabilizations were used in the 
fictitious domain method of [12]. However, the set of faces that are stabilized are slightly different 



here and we choose s in equation (2.34) in different ways depending on if the interface cuts the 
domain boundary or not. 



3. Analysis 



In this section we will show that the finite element method presented in Section [2~4| has optimal 
convergence order. We begin by proving the following consistency relation for the finite clement 
formulation (2.27). 
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Lemma 3.1. Let (u,p) £ WxV be the solution to the boundary value problem (2.1) and (uh,Ph) 
be the solution of the finite element formulation (2.27). Then 



A h (u - u h ,p-p h ;v h ,q h ) = e u J(u h ,v h ) + e p J{p hl q h ), V(v h ,q h ) eW/,X Vh- 



(3.1) 



Proof. First, note that since k\ + K2 = 1 we have 

[abj = {a} [6] + [a] (6), (3.2) 

hence we can write 

j \(n ■ ( M e(«) - pl))v} ds = J^{n ■ (fie(u) - pi)} [«] ds + J \n ■ ( t ,e(u) - pi)} (v)ds. (3.3) 

Using the interface conditions for the normal stress and the shear stress, equation (2. Id) and (2.2), 
we have that 

J [n • (/ie(u) — pl)j (v)ds = J (jn(v-ri)ds. (3-4) 

Now, multiplying (2.1 ) by a test function (vh, qu) € Wh x Vh and integrating by parts, using ( |3.3[ ) 
and (3.4), the boundary conditions (2.1e) and (2.3), and that u is continuous (2.1c) we get 

(3.5) 



ah{u,v h ) + b h (u,q h ) - b h (v h ,p) = L h {v h ), 



and the claim follows. 



\Vh\\\h 



We introduce the following mesh dependent norms 

III^IH 2 = |||i 1/2 e(i;)||g, niU n a + ||{n • £ 1/2 e(«)}||* i/ 2lh ,r + ||M 1/2 HII Wr 
+ \\n ■ /i 1/2 e(«)ll-i/2 lfc ,*i + \\p}' 2 v\\l,2,h,an V« € W + W h , 
■■ \\\v h \\\ 2 + J(v h ,v h ) Vv h eW h , 

: lll^lll 2 + llM- 1/a ?HS, ni un a + U»r 1/2 ti}\\ 2 - W 
+ U- 1/2 qf-i/2, h ,m V(«,g) £ (W + W h ) x (V + Vfc), 
IMII 2 + llA* -1/2 «fcllo,n h , 1 un hia + J(lh,qh) V(v h ,Qh) € W/, x V h , 
where /2 1 / 2 = /ij/i} -1 / 2 , the norms on the trace of a function on T are defined by 

II 2 



[v, 



(Vh,qh)\\\h = 



ll«lli/a,fc,r = h K l W v \\l,rriK> 



IMI-i/W = S h «W v Wo 



TDK' 



□ 



(3.6) 
(3.7) 

(3.8) 
(3.9) 

(3.10) 
(3.11) 



KelCr 



and similarly for the trace of a function on d£l. Here /Cr is the set of elements that intersect the 
interface but when a part of T coincides with an element edge only one of the two elements sharing 
that edge belongs to /Cr- Note that 

(w,v) r < |Mli/2ArNI-i/2,fc,r, (3-12) 
{w,v) dn < \\w\\i/ 2 ,h,dn\\v\\-i/2,h,dn- (3.13) 

We will need the following inverse inequality when proving the inf-sup stability of the finite element 
method. 

Lemma 3.2. Assume that K is an element in /Cr such that, for i = 1 or 2, \K D = ai_Kh 2 K , 
where oh k > and let \Tr\K\ = Jk^k- For any function v € Wh, the following inverse inequality 
holds 

„2 / «?7if n,/..M,2 (3 M) 



h K \\n t e(v) ■ n\\' z ornl? < 



Oii.K 



\ e ( v )\\o,Knni- 
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Lemma 3.2 follows using that the functions in Wh are linear, cf. |17j . 

We also state two trace inequalities that we need for proving an approximation result. 

Lemma 3.3. Let K £ K. and v G H (K). There exists constants C and C such that for s G 



K 



o,rni? 



< C(h 



K 



\0,K 



v\\l, K ), 

h K s \\v\\l dK < C(h-'- s \\v\\l K + h£'\\v\\*). 



(3.15) 



Under Assumption 1-3 the first trace inequality follows from Lemma 3 in |10j and a scaling argu- 
ment. The second trace inequality follows from a standard trace estimate, see [HI Theorem 1.6.6]. 
We will also need the following estimates: 



Lemma 3.4. Let K e K and v G Then, 



\Vv\\l K < Ch- K 2 \\v\\l K , 

\\v\\l m < Ch ~K\Hl K , 



\ v \\a.rr\K 



<CKg\\v\\l 



\0.K- 



(3.16) 
(3.17) 
(3.18) 



The inverse inequality (3.16) follows from |18l Lemma 4.5.3] and the trace inequalities follow from 
Lemma 3.3 and the inverse inequality (3.16). 

3.1. Continuity 

We begin by showing the continuity of a/,(-, •) and then we prove the continuity of Ah{-, •; •, •). 
Lemma 3.5. Let u G W + Wh- There exist a constant C con t such that 

a h {u,v h ) <Ccont\\\u\\\\\\v h \\\ Vd/,6% (3.19) 



Proof. For each K G /Qr let |rni4T| = ^k^k and \K\ — axhk for some jk, ax > 0. The claim fol- 
lows by recalling the definition of a? l (-, •), applying the Cauchy-Schwarz inequality, inequality (3.12[) 



and (3.13) to the interface and boundary terms, and noting that {n-^e(u)} = {n-^ 1 / 2 e(it)}{/i} 1/2 



A| rn7? < (C r max \£ ) + D ) ^ A| M < (3.20) 

□ 

Lemma 3.6. Let u G W + Wh, Vh G Wh, P G V + Vh, q 6 V/,. There is a constant C con tA such 
that 

Ah(u,p;v h ,q h ) < C contA \\\(u,p)\\\ |||(i>h,g h )lllfc (3.21) 



Proof. 



Ah(u,p;v h ,q h ) = a h {u,v h ) + b h (u,q h ) + b h (v h ,-p) (3.22) 
= 1 + 11 + 111 (3.23) 



Term I. Using the continuity property of <%(•, •) (Lemma 3.5), we have 



I< C cont \\\u\\\ \\\v h \\\. (3.24) 
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Term II. Using the Cauchy-Schwarz inequality and the trace inequalities in Lemma [3^4] to bound 
the contributions from the interface and the boundary we get 

II < ||m 1/2 V • u|| ,fi 1 un 2 ||M" 1/2 5/ t ||o,n 1 ua 2 + ||M 1/2 M ||i/2,fc,r||M" 1/a {8/,}||-i/2,h,r 
+ \\v 1/2u \\i/2,h,dn\\v~ 1/2 qh\\ -i/2,h,en 

< C\\\u\\\ (ll/i-^^Ho.n.ua, + ||M~ 1/2 {«h}||-i/2,h,r + II/^VH-iamu) 
<C\M\^~ y2 Qh\Wa K1 ^- ( 3 - 25 ) 
Term III. 

1/2 



III < 



CIIKIII (|I^ 1/2 pIIo A uo 2 + ll{M}- 1/2 WII-i/2 ) / l ,r + U- X 'M-i/2,h,on) • (3-26) 



Summing the estimates of terms I — III and using the definition of the norms |||(-, -)lll an d |||('>')llk 
yields the claim. □ 



3.2. Inf-sup stability 



We will show that the finite element formulation (2.27) is inf-sup stable. Namely. 



Theorem 3.1. Let (v,h,Ph) G W/j x Vh- There is a constant C such that 

A h (u h ,p h ;v h ,q h ) + e u J(u h ,v h ) + e p J(p h ,q h ) 

sup rr-r- \Tn >C s \\\(u h ,p h )\\\ h (3.27) 

(v h ,q h )eW h xV h \\\\ v h,qh)\\\h 

First, we need to show the coercivity of dh(-, •). 

Lemma 3.7. There exists a positive constant C coer such that 

Ccoer-IIKH^ < a h {v h ,v h ) + J(v h ,v h ) Vv h € W h . (3.28) 

Proof. Let Ki = K Dili- By the definition of a^(-, •) we have 

2 

a h (v h ,v h ) = W^^VhMln, + ( A rlrn7?ll M llo.mx ~ 2 ^ n ' ^Ofc)}, HDrnA 

i=1 A'e/Cr 
2 

+ Y Y {*n\d£mKi\\ v h,i\\o,dnnKt - 2 ( n 1 ^(Vh,i),v hl i) a n nKi ) ■ (3.29) 

i=l KGKan 

For each if 6 /Cr let |rnif| = jk^k, \K\ — axhk, and = a^xhk- Using the Cauchy-Schwarz 
inequality and the geometric-arithmetic inequality we have 



Y ( A Hrnidl M IIojyia _ 2 ({ n ' M 6 ^)}' M)rnK 

Ae/Cr 



> 



Yl (^rlrng- fe^H lv h j |lo.mi?~Z! ^ ^ IK" • M^K,*)!! 2 , rng - 

\ 7A / i 

A'e/Cr 1=1 AS/Cr 

(3.30) 

Here, we take 

(l + fl^/thl 2 A (l+i?)^ 2 ^72 2 

5l ^ = AcV, ' = Ac7 2 ' (3 ' 31) 
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with oii = di^K and 7,; = jk when T n K intersects the element boundary exactly twice, and 
ati = a F , 0L2 = ax, 7i = 7t ; 72 = Ik when TDK coincides with an edge shared by element 
T C Oi and K C ^2- With our choices of «i, K2, and Ar, using the inverse inequality in Lemma|3]2j 
Hi = /2| /2 {m} 1/2 , and that > C we obtain 

2 

a/,(«h,o h ) > 5^(1 - A)\\^ /2 e{v h M. n , + ^C||{n • M 1/2e ( u ft)}ll-i/2,fo,r 

2 

+ £)||{/i} 1/2 ||i /2A r + X] H ( A n|annK < ||«/i,i||o,8nnic j ~ 2 ( n ' Vi*{vh,i), v h,i)dar\K % ) ■ 

i=l KelCau 

(3.32) 

Applying the Cauchy-Schwarz inequality, the geometric-arithmetic inequality, and the inverse 



inequality (3.34) on the boundary terms we get 
2 



2 



an \ 1 ian,K 



2 



where Jdn^hic = \dSl n if|. For each K £ /Can the following inverse inequality holds 

h K \\n ■ €(w M )||j» ianrts . < T^lle^fc.OHg.x. (3.34) 

For elements if g /Can such that if n fl; 7^ 0, if <f_ Hi, let J-'k k* be the set of all faces that has 
to be crossed to pass from K to K % C and iVp the number of such faces. Our assumptions 
guarantee that such an element K l exists and that there are a bounded number of faces in J- F ,K i - 
We use the same idea as in [11] and write that 

e(vh,i)\K = e(vh,i)\K* + tf[wF- e(v hii )}pn F , (3.35) 

where 5 — ±1 with the sign depending on the orientation of n F so that the equality holds. We 
then have 

||e(« M )||§ )K < 2 I ^Mv h MIk*+Nf^ E II I"^«m)Uo,f ) , (3.36) 

where we have used the Cauchy-Schwarz inequality and the geometric-arithmetic inequality. Due 
to Assumption 1 (quasi- uniformity) we have that there exists a constant c q = max fppji jpjfc 
and hence 

\\e(v hti )\\lK < 2c 9 I HeK^Ho,^ + N F h £ || \n F ■ e(v Kl )j F \\l F I . (3.37) 

Let N K i be the number of elements K € /Can, K $L Vti that have if 1 as the closest element 
completely in fij. Taking 

2c,(l + F)M»7gn,K . . . ,„ 
Odn,i,K = ^ « = 1,2, (3.38) 

iiQiK 
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when K (IT intersects the element boundary exactly twice and otherwise 

,/2 



Son,-, 



(1 + F)wl nK 



i = 1,2, 



(3.39) 



we have from our choice of Agn (equation (3.33 1) the inverse inequality (3.34), equation (3.37), 
and ( "p™-" ) > C that 



(v h , v h ) >^2(1-A- E(l + max(A^O))lk 1/2 ^M)llo,n, + BC\\{n ■ ^ 2 e(v h )}\\\ 



D\m 1/2 W || 1/2A r + FC\\n ■ ^ 1/2 e{v h )\\ 2 _ l/2 ^ 



on 



+ Ghi 1 l 2 v h \\\ /2hm ~Y,EN F J2 J2 IHh\\ln F -e(v hti ) 



If No,f- 



Note that for each face in F e Tk,k* we can write 

=e(«fe,i)|jf- +ln F ■ e(v h:i )} F n F , 
V(v?j,i)| K + = V^/j.i)!^- + [nj? • V(v, l ,i)] F n F , 

and that each of the terms in equation ( 3.41[ ) are constant 2x2 matrices. Hence, 



J] \\ln F -e(v lhi n F \\l F < Wl n F-V(v h ,i)l F 



2 

0,F- 



Finally, since EN F < j 



(1-A)N F 
+max(jV K i ) 



< 1 we have 



/2,h,r 



(3.40) 



(3.41) 



(3.42) 



J2EN F J2 lM h \\in F -e(v Kl )l F \\l F <J(v hl v h ) 



KGKgn, F£F K K i 

Knfii^ti. 
KgLtii 



and the claim follows. 



We will need the following technical lemma. 
Lemma 3.8. Let ph — (ph,i,Ph,2) G V^. There is a constant C such that 

Wth 1/2 PhA\l,n h , z < C (\\lii 1/2 Ph,i\\o,u, h , { + J(Ph,Phj) 



(3.43) 



□ 



(3.44) 



Proof. For elements K £ K-i that are not entirely in f2j, let Fk.k 1 be the set of all faces that has 
to be crossed to pass from K to the closest element K l c Sli and N F the number of such faces. 
Assumption 4 guarantees that such an element K % exists and since the mesh is assumed to be 
shape regular there are a bounded number of faces in Tk.k' 1 ■ We can write that 



Vh 



v,i\K =Ph,i\K* + X! s l n F ■Vph,il F n F ■ (x - a F ), 



(3.45) 



FeF R 



where S = ±1 with the sign depending on the orientation of n F so that the equality holds and 



a F is the center of gravity of F. Taking the square on both sides of identity (3.45), using the 
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Cauchy-Schwarz inequality and the geometric-arithmetic inequality, we get 



lbh,i|lo,iC 



<2 [ ^\\Ph,i\\o,K* + N F E ]^fllI^- V PMlF^-(x-a F )||^ F 



K.K % 



< 2c 9 \\phX t K* + N F E h *\\ ■ VP^F \\l,F > ( 3 - 46 ) 



Fer T< 



where the last inequality follows due to Assumption 1 (quasi-uniformity). Let be the number 
of elements in /Cp that have K l as the closest element completely in f2j. Summing over all elements 
K G K-i and using equation ( |3.46[ ) for element that are not entirely in £li we obtain 



Ik 1/ V 1 <llo,n fc . < < E 0- + 2c i N ^) Ik 1,2 PhAl,K> 

+ 2c 9 max(iV F ) £ £ ^ h 3 \\ [n F ■ Vp M ] F ||o> (3.47) 

KeK^KgLQi Fer K K i 

□ 

To prove the inf-sup stability of •) we use some of the ideas in [19] . Introduce the piecewise 
constant function 

P = \ Ml '^d"l-i ""o ( 3 - 48 ) 
[ on fi 2 - 

Let M =span{p}. For any p^ € we can write 

Ph = Po +Pft,o> Po € Mo, P^o e M^ . (3.49) 
The functions in M^ satisfy (p^ , = 0, i = 1, 2, see [19]. 
Lemma 3.9. For any po € Mo, i/iere exists Vh,o € suc/i i/iai 

6fc(«ft,o,Fo) > Ci : p ||^ _1/2 po|lo^ 1 uo 2 ' IlkMlIU < C 2 , Po ||^ -1/2 po||o,niun 2 - ( 3 - 50 ) 

Proof. Let po = po, then (po, l)niun 2 = 0- Let I (pa) be the continuous piecewise linear 
approximation of po which differs from po only in elements K S -Kr- Let a be the constant for 
which (I (pa) — a, 1)q — and let qn = I(po) — ol. Since the underlying finite element spaces are 
inf-sup stable there exist Vh,o & n C(fii U il 2 ) with u^olan = such that 

b h (v h fi,qh) r ,, || ,„ 



We have 



II Vf ft,o||o : Oiuo 
b h (v hfi ,Po) b h (v hfi ,q h ) b h {v hfi ,Pa - qh) 



II Vvft i o||o,n 1 ua 2 ||Vuft,o||o,niun 2 II Vv^ollo.niuOa 

> C||<3^||o,niu^ 2 _ ^2||Po - ?/i||o,fiiun 2 

> C , ||p ||o,n 1 un 2 - (V^ + C)IIp"o - Shllo.niuns 

> 6 - (6 + V2) — r— -j: ||Po||o,fiiun 2 

V ||Pollo,fiiun 2 / 

> (c-ch 1 ' 2 ) Upollo.nxu^, (3-52) 



where we in the last step have used that 



\{I(p ) - po, l)niun 2 | < WHPo) - Polio, Qiun 2 , g ^ 

Pllo.fJiuna ||l||o,fiiun 2 
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and hence 

||po - gfe||o,niun 2 < q II-^Po) ~ Pollo,niun 2 < ^ti/a 

IIpoII 

From the definition of Mq one can see that 



(3.54) 



II A* 1/2 Po\\o,n 1 un 2 = c (^ty\\Po\\l,a 1 un 2 , b h (v h ,o,Po) = C(iJ,,Q)b h (v hfi ,Po), (3.55) 

with C(/i,fl) = Ml |ni|-i+fnl"-i — ■ ^ e can choose Vh,o so that equation (3.52 1 is satisfied and 
HVuft.ollo.niufia = ||po||o,niun a and obtain 

bh(vh,o,Po) = C(fj,,n)b h (v hfi ,p Q ) > cC(fi, fi)||po||o,niun a = c llM _1/2 Po|lo,niun a - 
We have 



C{n, Q) > min 



\Oi 



V|fii| + |n 2 



r 



(3.56) 
(3.57) 



and 



HK,olll ^ C^^cl|Vw/,,ollo,n 1 un a = C , /iJ n / a 2 x ||p |!o,n 1 uo 2 < CC 1/2 ||^ 1/2 Po||o,o 1 un 2 - (3.58) 

Finally, note that for Vh,o = { v h,i, v ha) e Wh H C(Qi Uf^) we can choose u^.i = v/^j in fi/^nfij, 
j ^ i so that Hlw^ollU < C|||«ft,olll- □ 



Lemma 3.10. For any pj[ a — {Phi,Phi) G -^To ^ ere exists Vh & Wh such that 
and \\\v h \\\ h < C||^~ 1/2 p^||o,n h l uf2 h , 2 - 



(3.59) 



Proof. Let a be the constant for which 

iPhfi - l )oj Ki = 



(3.60) 



and let qh,i = p^ j — a. Using the fact that the underlying finite element spaces are inf-sup stable 
with a uniform constant on any polygon of shape regular elements which has no clement with 
two edges on the boundary, see Brezzi-Fortin [2D], Proposition 6.1, Page 252, we can for each 
qh = (qh,i,Qh,2) choose Vh,n t = Ki,»ft,j) € Wh with v hA € Wh,i such that v h ^ = on fth,i\^h,i 
and on dfl and Vh,j — for j ^ i so that 



bh( v h,niiQh) . ^,1, -1/2 



Using the inverse inequality (3.17) and that supp(w^ i si i ) = Wh,i 
which together with Korn's inequality [21] Eq. (1.19)] yields 



III*/., 



UK* 



+ J(vh,ni,v h ^ t ) < C\\\v h ,Qi 



(3.61) 



(3.62) 



(3.63) 



We can choose v/j^i so that equation (3.611 and (3.63) are satisfied and |||i'/i,n i \\\h = ||/^ ' 9/i,i||o,n h ,, 
We then have 

bh(vh,Qi,qh) > C q \\^ 1/2 q h A\\o,n h J\N 1/2< lh : i\\o^ h]1 - (3.64) 
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Lemma |3.8| then yields 

ll/^V'lliW < C (M 1/2 q h ,i\\l Uhli + ■%,,«,)) < C (0^^^) + J(q h ,q h )) 
Note that b h (v h ^.,q h ) = b h (v h ^ v p^ fi ), J{q h ,qh) = J(Ph,o>P~h,o) and 

\(f ii 1/2 a,l) Uhi \ _ \{Lh 1/2 Ph^ l )nA^ h ,,\ \\^ 1/2 Ph,i\\o,n h ,, II l||o,n,W., 

\\ni u > t llillL,,. : PIIL^ 



ImT 1/2 -l = 



(3.65) 



(3.66) 



Using equation (3.66) we get 



-1/2 i|2 ^ II -1/2 _L ||2 || -1/2 ||2 

-V2 T1 X.||2_ | j 



> IK V Xillo,n hii 



ll 1 llo,n M ll 1 llo,niW,< 



Il0,w h ,i 



We assume ||l||o,n<\w*, 

IIa'- 



From equation (3.66) we 
■w/i,n 2 we have 



.., _ . = ch 1 / 2 and hence equation ( |3.65 ) and (3.67) yield 

1/a PMllo,n M < C(l - ^ 1/2 ) _1 (C^K^o) + J(Pm,Pm)) ■ 

also obtain Hlw^nJIU ^ Cll^^PiJo,^,, • Finally, taking u h 



2 



and ||K||| h < CELi II ^ i^illo,n M - 



(3.67) 

(3.68) 

w/i,ni + 

(3.69) 
□ 



Lemma 3.11. For any ph € i/iere exists Vh € and constants C%, C3 > 0, and C2 > suc/i 
f/iaf 

b h (v h , Ph ) > Ci||/x~ 1/2 j>h||o,n )k !un h a - C 2 J{p h ,Ph), \\\v h \\\ h < C s \\n~ 1/2 Ph\\ofl h ,iun hJ ,. (3.70) 



Proof. If p h G Vh is piecewise constant, i.e. ph — Po where po € M , the proof follows from 
Lemma 3.9 with C2 = 0. Otherwise, we have ph = po +p^ , where po £ Mo and pj^ G Mq h . Let 
Vh be such that Lemma 13.91 is satisfied and Vh such that Lemma 13.101 is satisfied. For a > 0, 



define Vh = Vh,o + ctVh- Note that bh{vh,Po) = 0, since Vh vanishes on T U dVL and po is constant 
on each subdomain i = 1,2, and 

\bh(v h) o,pj; fi )\ = |(V-i; M ,p^o)n 1 un 2 | < C|||7;, li o|||||^ 1/ W,ollo,n 1 un 2 , (3.71) 

since Vh,o is continuous and vanishes on dtl. Also, J{ph,Ph) = J(Pho> 

bh{vh,Ph) = b h (v h .Q,p Q ) + b h (v h ,o,Phfl) + ub h (v h ,po) + ab h (vh,Ph,o) 

> C r i,pollM _1/2 PQ|lo,fiiun 2 ~ c ' c 2,p ||M~ 1/2 Po||Q,n 1 ur2 2 ||M _:1 ' /2 Pft ) ollQ,^,iun il , 2 
+ a ( C hPiJ^ 1/2 Phfi\\o,n h ^un h , 2 - C2,p^ J(Ph,o,Ph,o)) 

> Ci||/i- 1 / 2 po||o,n hll unH >a + ( aC h P l» ~ c a) llM~ 1/2 Ph,ollo,n fc , 1 un hia 

- uC 2Ko J(pl ,pl ) > Ci||M- 1/a Ph|lS,n hil un fc , a - C * J (Ph>Ph)- (3.72) 
For sufficiently large a. Finally, we also have 

Hb/illU < |||«mIIU + a lll*fc||U < C\\v~ X/2 Ph\\o,Q. h iun h 2 - (3.73) 
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□ 



We are now ready to prove the inf-sup theorem. 
Proof, (of Theorem 3.1) Note that if ph € Vh is constant, ph — since (fi~ 1 pi l , l)f2iun 2 = 0. 
Letting Vh = Uh and using the coercivity of a^(-, •) we have 

A h (u h ,p h ;u h ,Ph) + e u J{u h ,u h ) +e p J(p h ,p h ) = a h {u hl u h ) + e u J(u h ,u h ) 

min(l, e u )\\\(uh,Ph)\\\h> (3.74) 



and hence the proof follows. Otherwise, let Vh be such that Lemma 3.11 is satisfied and 

and 



Then, using the coercivity and continuity of a/j(-, •) (Lemma 



3.5 



\W 1,2 Ph\\a,n h .iun h ,2- 

Cauchy-Schwarz inequality, and the stability of bh(-, •) (Lemma 3.11 ) we have for a > 
A h (u h ,p h ;u h - av h ,Ph) + e u J{u h ,u h - av h ) + e p J(p h ,p h ) = 

a h {u h ,u h ) + e u J(u h ,u h ) - a (a h (u h ,v h ) + e u J(u h ,v h )) + ab h (v h ,p h ) + e p J(p h ,p h ) > 



3.7), 



C coer min(l,£„)|||'U/ l |||£ - amax(C cont , e u ){\\\u h \\\ + {J{u h ,u h )) 1/2 )(\\\v h \ 

+ a (c 1 \\ f j,- 1 ^p h \\l Qh lUQh2 - C 2 J(p h ,p h )^ +£pJ(p h ,Ph) > D\\\(u h ,p h )\\\l. 

with D = min(Di, D 2 , -D3), where 

Di = (Ccoer min(l, e u ) - amax(C cont , e u )/5) > 0, 
D 2 = a (Ci - max(C co „ t ,e„)(5) > 0, 
D 3 = (e p - aC 2 ) > 0, 

provided 5 and a are sufficiently small. Finally, the proof follows using that 

\\\(u h - av h ,p h )\\\ h < \\\(u h ,p h )\\\ h + a\\\v h \\\ h < (1 + a)\\\(u h ,p h )\\\ h 



(J(v h ,v h )) 1/2 ) 

(3.75) 



in equation (3.75). 



(3.76) 

(3.77) 

□ 



3.3. Approximation properties 

In this Section we will show that the spaces Vh and Wh have optimal approximation properties 
on i? 1 (i7i U 5^2) and [H 2 (Vli U f^)] 2 , respectively, in the energy norm. In order to construct an 
interpolation operator we recall that there is an extension operator £* : H s (rii) —> H s (il), i = 1,2, 
s > 0, such that £^Wi\n i = Wi and 

||£>ilkfi < C|K||,,n 4 Vu* e i?*(fii). (3-78) 

See [55] for further details. Let tt^ : H s (il) —> V^ be the standard Scott-Zhang interpolation 
operator [33] and recall the stability property 

h h w\\r,n < C||HU,n, 0<r<min(l,s), \fw£H s (fl) (3.79) 

and the approximation property of the interpolation operator 

\\w - TT h w\\ r . K < Ch S x r \w\ sM(K) , 0<r<s<2, \JK e JC, Vw G H s (n), (3.80) 

where AT(K) is the set of elements in JC sharing at least one vertex with K. We define 

nliWi =ir h £fwi\„ hti VtOi € H'(SU) (3-81) 

and for w = (wi, w 2 ) with U»j|n 4 € H s {p.i) we define 

Kw = «,i^i,<,2^2)- (3.82) 

We will use the same interpolant for velocities and pressures. For the velocities s = 2, Vh.o — Wh,o, 
and TT* h ■ : H 2 (Q,i) -> Wfc,i, i = 1, 2 while for the pressures s = 1, V/,,o = and Tr h i : H 1 ^) 
Vfc,i, i = 1,2. In the norm |||(-, -)|||, we have the following interpolation error estimate: 
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Lemma 3.12. It holds that 

\\\{v--K* h v,p-KtpW < h^CJ^Mln^+Cpll^Wi^unJ- (3-83) 



Proof. Recall the definition of the norm |||(-,-)lll (equation (3.8)). The interface and boundary 
contributions can be estimated in terms of element contributions by applying the trace inequal- 
ity in Lemma |3.3| Then, for the element contributions, applying the approximation property of 



the interpolation operator (3.80), and finally using the stability of the extension operator equa- 



tion (3.78) yields the desired estimate. rj 



3.4- A priori error estimates 

Theorem 3.2. The following error estimate holds 

|||(«-u h ,p-p h )lll < Ch (llM^lkfixun, + ||/i" 1/2 p||i,n lU n 2 ) (3.84) 

Proof. We have 

\\\(u-u h ,p-p h )\\\ < \\\(u-Trlu,p-TTlp)\\\ + |||««- u h ,w* h p-p h )\\\ h . (3.85) 
Here the first term can be estimated directly using the interpolation error estimate (Lemma 3. 12| ) 



(u - n* h u,p- n* h p)\\\ < Ch {\\^u\\ 2 , ni un 2 + ||M~ 1/2 Plli,n 1 un a ) • (3-86) 



Turning to the second term we use the inf-sup condition (Theorem 3.1 ) followed by the consistency 



relation, Lemma 3.1 to get 



\\\(% h u-Uh,Tr h p-ph)\\\h 

, n -i A h (ir* h u - u h ,Ti* h p-p h ;v h ,q h ) + e u J{ir* h u - u h ,v h ) + e p J(-K* h p - p h ,q h ) 

S ^ s SU P nTT \TTi 

(v,g)eW h xV h II \K v hi lb.) 1 1| h 

^ r>-\ A h (n* h u-u,Trlp-p;v h ,q h ) + e u J(Tr* h u,v h )+e p J{iT* h p,q h ) 

S^ s Slip 777- rjji (d.SfJ 

(v,q)eW h xV h \\\{Vh,Qh)\\\h 
Using the continuity of Ah{-, •; •, •) and that from the Cauchy-Schwarz inequality we have 

J(n* h u,v h ) < J(7r* h u,7r* h u)^ 2 J(v h ,v h )^ 2 < J^u.^uJ^IHK.gOIIU (3-88) 
and similarly for J(Tr^p,qh) it follows that 
Uh,n* h p - Ph)\\\h < 

C; 1 (c contA \\\(u-ir* h u : p-ir* h p)\\\ + e^tu^u) 1 / 2 + e p J(n* h p,7T* h p)^ 2 ) . (3 . 89) 

The first term is estimated using the interpolation error estimate. For m G H 2 (tti) and £ 2 u = 
(£ 2 iti, E2U2) we have that 

J(tt^u 7 tt^u) = J(£ 2 u — tt^u, £ 2 u — 7t£m) 



i=i Fejr r ; 
2 

<Y, C ^Y, (^? u * - <M\Ik + h 2 \\£ 2 Ul - km\Ik) 

i=l KeK 

2 2 

< Ch 2 J2»iJ2 W £ ? u iWlK < C^E^H^II^, (3-90) 



»=i fce/c i=i 
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where we have used the Cauchy-Schwarz inequality, the trace inequality in Lemma |3.3| the ap- 



proximation property of the interpolation operator (equation (3.80)), and finally the stability of 



the extension operator (equation (3.78)) 



Third term in equation (3.89) can be estimated using the inverse estimate 



h\\ [n ■ Vtt*^] \\l F < C\\V7rl iPi \\ 2 0K+UK -, 



(3.91) 



where Kp and K F are the elements sharing face F. The inverse estimate together with the stability 



property of tth equation (3.79) and the stability of the extension operator (equation (3.78)) yield 



J«p,<p)=£ C ^\\ Inp-^liPi] \\o,f 

2 2 

< ch 2 \h Y E \\^si Pi \\i KtUK - < ch 2 J2 ^ E \whi\\l K 

i=l FgJr.i ' *=1 kelC 

2 

<Ch 2 J2^\\ Pi \\l ni . (3.92) 



Collecting the estimates ( |3.85| ), ( |3.86[ ) ( |3.89[ ), ( |3.90[ ), and ( |3.92| ) the theorem follows. n 
An L 2 -estimate for the velocity can be proven assuming additional regularity and using the Aubin- 
Nitsche duality argument, following the proof of O Proposition 11]. 



4. Estimate of the condition number 



Let {^Pi}^ and {Xi}£i be a standard finite element basis in Wh and Vh, respectively. Let A 
be the stiffness matrix associated with the formulation ( 2.27 ). Matrix A has dimension (Ni +N2) x 
(Ni + N 2 ). For the Euclidian norm of a vector X £ R N we use the notation |A|^, = J2i=i ^H- 
We recall that the spectral condition number k(A) is defined by 



n{A) = \A\ N \A- X \ N . 



(4.1) 



Here N — (Ni + N 2 ) and \A\ N = sup| X[jv=1 \AX\ N for A G R NxN . The expansion it/, = 

Si=i ^iVi an d P/i = X)i=i PiXi define isomporphisms that map tt/j e W/i to U G K^ 1 and 
Ph £ V?,, to P £ , respectively. We have for V € 1* being the concatenation of U and P the 
following estimate 

c\h~ x (IMIo.n^un^ + H^lkn^un^) < |V|jv < c 2 h~ 1 (\\Ph\\o,n hA un h:2 + \\uh\\o,n hA un h , 2 ) • 

(4-2) 

To derive an estimate of the condition number we first prove a Poincare type inequality in 
Lemma 14.11 and an inverse estimate in Lemma 14.21 Then the condition number estimates follows 
from these two lemmas and the approach in |23] , 



Lemma 4.1. The following estimate holds 

max ) A*min 

for all (v h ,q h ) £ W h x V,,.. 



,1/2 



l/2\ 



\{vh,Qh) 



(4.3) 



Proof. We have that ||g^||o.o h ! 

r,r 1/2 \ 



< (J-ULWfJ. 1/2 qh\\osi h , 1 un h . 2 - We need to show that ||u h ||o,n hil L 



j^ftlllh. Let <fi be the solution of the dual problem 

- V • e(4>) = v h inn, V • <f> = in Q, (f> = on dfl. 



(4.4) 
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We assume that we have elliptic regularity so that 

IMl!,n<C n |K||e> in . (4-5) 

Multipliyng the dual problem with Vh, integrating by parts, and using the Cauchy-Schwarz in- 
equality we get 

IKIIo,niun 2 = ( e ( t, / J ) 5 £ (^))o,o 1 un 2 - (M,n-e(0)) r - (v h , n ■ e{4>)) m < 

C (||e(wfc)||o,n 1 un a + l|[»h]|li/a,fc,r + \\ v h\\l/2,h,an) (ll 6 (^)llo,n 1 un a + IK <<l>)\\-i/2,h,r) 

(4.6) 

Using that \\n ■ e(4>)\\o r < Cr||<^»||| Q and inequality ( |4.5| ) we have 

<c|Kllo (4.7) 

Thus, 

(l|eK)||^ ni un 2 + \\[vh]\\l /2 , h ,T + Klli/^an) < ^/vlfilKIII- (4- 8 ) 
Following the proof of Lemma |3.8| we can show that 

\Mo,n h<i < C ||^||o^uo 2 + E E ^III^'^mIfIIo.f I (4-9) 



i=i 

Finally, we have that 

2 



V 



i=l Fefr,; 



E K,<llo,n M < O (/ViflKIII + ^L\j{v h ,v h ))^) < C^ 2 \\\v h \\\ h . (4.10) 



i=i 



Recalling the definition of the norm |||(-, -)\\\h, the claim follows. rj 

Lemma 4.2. TTie following estimate holds 

\\\(v h ,Qh)\\\h < C^ 1 max(^ / a 2 x ,/i m [ n /2 ) {hhWo^un^ + IKIkn^un,^) (4.11) 
for all (v h ,q h ) € W ft x V h . 

Proof. Note that \\n~ 1/2 qh\\o,n hA un h , 2 < ^n\\Qh\\o,n hl ua K , a - Using the Cauchy-Schwarz in- 



equality and the trace inequality (3.17) we have 



2 2 

J^,*) <E^ E Ch 3 \\tVq h 4 F \\l F < ^Ch*Y, E H V( 7mIIo,k 

i=l F6.Fr,; i=l ife/Cj 

2 

< /Vin^E E \\9h,i\\U < ^rlnC\\q h \\l nh lUnh2 . (4.12) 

In the same way we obtain 

J(v h ,v h ) < Hn^Ch s - 3 \\v h \\l nh iUnK2 . (4.13) 
The standard inverse inequality (|3.16[) yields |l/i 1/2 e(«/ I )|lo,n 1 un 2 < CV^Vm^ll^fclkniuna- The 



Lemma follows using the trace inequalities (3.17) and (3.18) on each of the interface and boundary 



contributions to □ 
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Theorem 4.1. The following estimate of the spectral condition number of the stiffness matrix 
holds 

K (i)<Cma%L^ m i n )ft" 2 . (4.14) 

Proof. We need to estimate |^4|jv and |^4 _1 |at. Let V € R N and W € be the vectors containing 
the coefficients corresponding to (ufc, (ft) € W/, x Vh and (uih,rh) & Wh X Vfe, respectively. Starting 
with |,A|jv we have 

. .... cav; wp^ 

|^V|jv= sup — j— | 

A h (v h ,q h ]w h ,r h ) + e u J{v hl w h ) +e p J(q h ,r h ) (a-\k\ 

— SUp . (4--L0J 

(tu h ,r ft )eW h xV h l^liV 

We now use the continuity of A/j(-, •; •, •) together with that \\\(vh, Qh)\\\ < C*\\\{ v h,qh)\\\h and the 
Cauchy-Schwarz inequality to obtain 

A h (v h ,q h ;w h ,r h ) + e u J(v h , w h ) +s p J(q h ,r h ) < 

Cmax(C*C cont A,e u ,e p )\\\(v h ,q h )\\\h\\\{wh,rh)\\\h- ( 4 -16) 



Lemma 4.2 and equation (4.2 1 then yield 



A h (v hl q h ;w hl r h ) + e u J{v h , w h ) +e p J(q hl r h ) < 

Cmax(/i max ,^~?- n )/i _2 (||?/ l ||o,n h , 1 un fc , 2 + \\vh\\o,n hA un hil ) (\\rh\\o,n hA un hi2 + \\iDh\\o,n hA un K2 ) 
< Cmax( Mmax , M -L)l^|iv|Wk. (4.17) 



Thus, we have the estimate 



\A\ N = sup N < Cmax(/x max , fij a ). (4.18) 
vm N \ v \n 



Next we turn to the estimate of \A x \n- Using equation (4.2 ), Lemma 4.1 and the inf-sup stability 



(Theorem 3.1| we get 

|V|jv < Ch^ 1 (\\qh\\o,n hA un h , 2 + ||'W/ l ||oj2 h . 1 un h . 2 ) 

< C max(^^ , ^i 2 )h- 1 1 1 1 ( v h , q h ) \ \ \ h 

<r™^(,Mi „- 1 /2^^-i m _ A h (v h ,q h ; w h , r h ) + e u J(v h , w h ) + e p J{q h ,r h ) 

^ ° max lMmax) Mmin j' 1 SU P [77777, 777771 

< Cmax(/i m ' m ,/jJ )/i \AV\ N ^ ^ = . (4.19) 

We conclude that \V\n < Cmax(/z max , /i~?" n )/i~ 2 | AV^at. Setting F = A~ 1 W we obtain 

< Cmaxt^^-Jjr 2 . (4.20) 



Combining estimates (4.18) and (4.20) of \A\jy and \A 1 1 jv the theorem follows. rj 
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Figure 3: Spectral condition number and the error in the pressure versus the stabilization constant e p for the 
continuous problem when h x = 0.0125. Left panel: The estimated spectral condition number versus e p . Right 
panel: The error in the pressure measured in the L 2 norm versus e p . The dotted line indicates the value of e p that 
we have used in the computations. 



5. Numerical examples 

We have shown that the proposed finite element method is of optimal convergence order and 
results in a well-conditioned equation system. In this section we present results using the proposed 



method (see Section 2.4) for numerical experiments in two space dimensions. We study the 
convergence rate of the numerical solution and the condition number of the system matrix for 
three examples. A direct solver is used to solve the linear systems. 

Unless stated otherwise, we report the size of the velocity mesh h x . The pressure mesh is twice 



as course. The parameters n\ and K2 are chosen according to expression (2.15) and the penalty 



parameter Ar is chosen locally according to expression (2.16). Dirichlet conditions for the velocity 



are imposed weakly and the penalty parameter Xqq enforcing the boundary conditions is chosen 



according to equation (2.22) 



Both the stabilization terms J{ph,<lh) and J(uh,Vh) are needed in order to have control of 
the condition number also in the cases when the interface is very close to a mesh line. The 
stabilization parameter e p = 1 and e u = 10~ 3 in all the examples. The errors are not sensitive to 
these parameters. Also, remember that we have defined fi = 2/x, in i = 1,2. 

5.1. Example 1: A continuous ■problem 

We consider a continuous problem presented in [3]. The computational domain is [0, 1] x [0, 1], 
the interface is a circle centered in (0.5,0.5) with radius 0.3 and fi = 2. The Dirichlet boundary 
conditions on <90 are chosen such that the exact solution is given by u — (20xy 3 , 5x 4 — 5y 4 ) and 
p = 60x 2 y - 20y 3 - 5. 



In this example we use a regular mesh. We choose Xqq according to equation (2.22) with G, 
F, and E such that Ago = Furthermore, we take — 2 and D — 0.05 in the expression for 
the penalty parameter Ap ■ The condition number of the system matrix and the error depends on 
these constants. However, we have not optimized these constants. 

In Fig. [3] we show the spectral condition number and the error in the pressure as a function 
of the stabilization parameter e p for h x — 0.0125. As seen in the figure the condition number 
of the system increases as e p decreases. Also, a too small e p results in a condition number that 
increases rapidly as the mesh size is reduced. However, for e p < 1, the error is not sensitive to the 
stabilization. Therefore, we have chosen e p — 1 in our computations. In this example the results 
for e u = 10 -3 and e u = coincide. Since the interface is not very close to any meshlines we have 
control of the condition number even when e u = 0. This is not the case in for example the last 
example in this section. 
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Figure 4: Convergence rate of the error in pressure and velocity for the continuous problem. Circles (o), Crosses 
(x), and stars (*) represent the new method, the standard continuous FEM, and the method in [3], respectively. 
Left panel: The error in the pressure measured in the L 2 norm versus the mesh size h x . The dotted line y = h x , 
shows the expected convergence order. Right panel: The error in the velocity measured in the L 2 norm versus the 
mesh size h x . The dotted line y = 2h x , indicates the optimal convergence order. 



The convergence for the velocity and the pressure in the L 2 norm is shown in Fig. [4] Since 
in this example neither the pressure nor the velocity have discontinuities we compare our method 
with the standard continuous finite element method. Compared to using standard continuous 
finite element methods we obtain just slightly larger errors for the pressure. However compared to 
the method in [3j (see Fig. 3 in [3] ) we obtain much smaller errors for the pressure. In Fig. [4] wo 
see the optimal second order convergence for the velocity in the L 2 norm but for the pressure we 
observe better convergence than the expected first order measured in the L 2 norm. In Fig. [H] we 
show the spectral condition number. The condition number using the proposed stabilized method 
grows as 0(h~ 2 ) just as it does for the standard finite element method, which is optimal. For 
a fixed mesh size the condition number of the system matrix using the proposed method is very 
close to the condition number of the system matrix using the standard continuous finite element 
method. We see that the condition number grows erratically with the mesh size when there is no 
stabilization for the pressure, i.e. e p = 0. However, we also see in the figure that the numerically 
estimated inf-sup constant in case e p = is essentially independent of the mesh size. Thus, our 
numerical results suggest that the inf-sup condition is satisfied in this case even when there is no 
stabilization. 

5.2. Example 2: Static drop 

Consider a circular interface T of radius R in equilibrium in the interior of a domain in two 
dimensions with fx = 2, a = 1 and vanishing on dfl. The exact solution is u = 0, p\ = 0, P2 = 
a / R. This corresponds to a circular fluid drop in equilibrium with the surrounding fluid. 

This problem has previously been studied for example in [53] using standard continuous finite 
element methods. In Figure[6]we compare the pressure approximation using the new method with 
the results obtained in [23] ■ In this case R = 0.5 and we prescribe the exact curvature k = 2. 

We use a regular mesh with h x — 0.025 in the velocity mesh. We choose Aon and Ar as in the 
previous example. From Table[T]we see that for the new method the magnitude of spurious currents 
and the error in the pressure are of the order of machine epsilon. However, using a sharp surface 
tension representation and standard continuous finite element methods the magnitude of spurious 
currents are large and may lead to unphysical movements of the interface. With standard globally 
continuous finite element methods the pressure either oscillates or is smeared out depending on if 
a sharp or regularized surface tension representation is used, see the two leftmost panels of Fig. [6| 
With the new method the discontinuous pressure is accurately represented even on coarse meshes. 
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Figure 5: Spectral condition number and the inf-sup constant. Circles (o), Crosses (x), and stars (*) represent 
the presented method with e p = 10~ 2 , the standard continuous FEM, and the unstabilizcd method (i.e. e p = 0), 
respectively. Left panel: Condition number versus mesh size. Right Panel: The inf-sup constant versus mesh size. 





\\u h \\oo 


Wp-PhWoc 


Condition number 


Regularized force 


o(io- 16 ) 


1.0129 


5.36 • 10 5 


Sharp force 


0.0126 


1.0164 


5.36 • 10 5 


New method 


0(nr 16 ) 


o(io- 16 ) 


6.74 • 10 5 



Table 1: Spurious velocities, error in the pressure approximation, and the spectral condition number for the static 
drop. Standard continuous finite elements with a regularized and a sharp approximation of the surface tension force 
are compared with the new method with e p = 10 _1 . 



Remark. We would like to emphasize that in order to get the magnitude of spurious currents and 
the error in the pressure of the order of machine epsilon even on course meshes it is important to 
use the bilinear form bh in equation (2.31|. Using bh in equation (2.30) results in larger spurious 
currents and errors in the pressure but the errors decrease with mesh refinement. Hence although 
the two forms bh in equation (2.30) and (2.31) are mathematically equivalent we obtain a perfect 



balance between the terms on the left and the right hand side of the variational form when bh in 



equation (2.31) is used 



5.3. Example 3: Couette flow 

We first consider a problem where there is a jump in the pressure due to different fluid vis- 
cosities. The computational domain is [0,L] x [—0.4,0.6]. The interface is the straight line y = 0. 
The viscosity 

a-( 200 y> °' (51) 

f = (3, 0) and the Dirichlet boundary conditions for the velocity are chosen such that the exact 
solution is given by 



u(x,y) = 
p(x,y) 



L — x 2 xy 
2L ' ~L 



2x 



Mi + A«2 



(5.2) 



We take L — 1. The interface intersects the boundary of the domain and max(N] K ) = 1. The 
penalty parameter Xgn is chosen according to equation (2.18) with E = 0.25, F = G = 0.005 and 
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Figure 6: Cross section of the pressure approximation for the static drop. The exact curvature k = 2 is prescribed. 
The dotted lines in all figures represent the exact pressure. Left panel: A Standard continuous finite element method 
and a regularized surface tension force is used. Middle panel: A Standard continuous finite element method is used 
with sharp representation of the surface tension force. Right panel: The pressure is approximated using the new 
finite element method. 



the constants in Ar are chosen as A — 0.3, B = D = 0.05. The condition number depends on 
these constants but we have not optimized these numbers. 

In Fig. [7] we show the approximation of the pressure using the new finite element method and 
the convergence of the error in the pressure and the velocity. We have as expected first order 
convergence for the velocity in the if^fiiUf^) norm and a bit better than first order convergence 
for the pressure in the L 2 (Ui U Q2) norm. 

Next we consider the problem presented in |2j. The computational domain is [0,L] x [0,H]. 
The interface is the straight line x = a. The jump condition |/iD(u) • n — pn\ ■ n = 1 is imposed 
at the interface. The Dirichlet boundary conditions on d£l for the velocity are chosen such that 
the exact solution is given by 

u(x, y) = (^iy( H - y), oV 

p(x,y) = -—x + X(x-a), (5.3) 

where X(x — a) = 1, if x > a and zero otherwise. We let L = 3, H = 1, a = 2, and // = 2 as in [2]. 

The approximation of the discontinuous pressure using the proposed method is shown in Fig. [8j 
The error in the velocity measured in the H 1 norm and the error in the pressure measured in the 
L norm are shown for different mesh sizes in Table [| The results in Table |] can be compared 
with the results in [31 Table 2]. The errors in the velocity are similar to the result in [2] but the 
errors in the pressure are much smaller using our method. 

In [9] we show the condition number as a function of the relative distance between the interface 
and the mesh line for different e u . We see that the stabilization for the velocity, J(uh,Vh), is 
needed in order to obtain a well conditioned system matrix independently of the location of the 
interface. 



6. Conclusions 

We have proposed a finite element method which offers a way to accurately solve the Stokes 
equations involving two immiscible fluids with different viscosities and surface tension. The in- 
terface that separates the two fluids can be represented either explicitly, for example as in the 
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Figure 7: Left Panel: Approximation of the discontinuous pressure in Example 3. The jump in the pressure is due 
to different fluid viscosities. The mesh does not coincide with the interface. There are 18 grid points along the 
x-axis in the pressure mesh. Right Panel: The error in the pressure measured in the L 2 norm and the error in 
velocity measured in the H 1 norm versus the mesh size h x . The dashed line represents y = h x . 




Figure 8: Approximation of the discontinuous pressure in Example 3 using the new method with e p = 1. There are 
35 grid points along the x-axis in the pressure mesh. The error in the pressure approximation for this simulation 
is shown in the third row of Table [2] 



h x 


\\u - Uh\\H 1 (n 1 un 2 ) 


lb-PfclU a (niun a ) 


Condition number 


1.88 • 1CT 1 
8.82 • 1CT 2 
4.41 • 1(T 2 
2.21 • lO- 2 

i.io- icr 2 


2.8742 ■ 1CT 2 
1.4137 ■ lO- 2 
7.6522 • 1CT 3 
3.6411 • 1CT 3 
1.8567- 10- 3 


6.8969 • 1(T 3 
1.6082 • 10~ 3 
4.6902 • 10~ 4 
1.0602 • 10~ 4 
2.7847 • 10~ 5 


1.55 • 10 4 
6.25 ■ 10 4 
2.14- 10 5 
9.25 • 10 6 
3.56- 10 6 



Table 2: Convergence study of the proposed method and the spectral condition number for Example 3. 
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Figure 9: The spectral condition number as a function of the relative distance between the interface and the mesh 



immersed boundary method, or implicitly as in the level set method. Our method allows for dis- 
continuities across the interface which can be located arbitrarily with respect to a fixed background 
mesh. 

We have used the inf-sup stable PI iso P2/P1 element and proven that our method is of 
optimal-order accuracy and that the stabilization terms J(uh,Vh) and J(ph,<lh) guarantee that 
the condition number of the system matrix is 0(h~ 2 ) independent of the interface location. We 
expect the method to be applicable also in three space dimensions. For higher-order elements, 
the stabilization terms J(uh,Vh) and J(ph,Qh) will include jumps of derivatives of higher orders, 
see [T5]. One can also include projection operators from [T3] into the stabilization to reduce the 
amount of stabilization and hence the constant in the error. The method we have presented is 
simple and robust and has properties that are very desirable, in particular for problems with 
moving interfaces. 
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